// Estimates the effect of natural gas prices on production facilities' emitting behavior and performs robustness checks
// Author Levi Marks levi.marks.1>at>gmail.com

// 0 // Preliminaries
global DATA "[Insert Path Here]"
set more off
graph drop _all

// 1 // Setup, trimming

use "$DATA/ghgrp_di", clear

// Merging prices (proprietary data - see ReadMe for instructions on how to incorporate)
merge m:1 facility_id year using "$DATA/Prices/ghgrp_prices"
drop if _m == 2
drop _m

// Dropping facilities whose DI production values (gas) differ substantially from their GHGRP production values (prod_mcf) by more than 20%
gen prod_ratio = prod_mcf/gas if year >= 2015
bysort facility_id: egen mean_prod_ratio = mean(prod_ratio)
gen trim_prod = 1 if mean_prod_ratio < .8 | mean_prod_ratio > 1/.8
drop prod_ratio mean_prod_ratio
drop if trim_prod == 1

// Only two remaining facilities in Alaska after trimming, which are effectively removed by fixed effects (one has six years after trimming, one has only one year)
drop if region == 1

// Dropping facilities with operations in either Colorado or ND
gen trim_co_nd = 1 if colorado != 0 | north_dakota != 0
drop if trim_co_nd == 1

// Trimming emission rate outliers at the facility level
drop if missing(rate)
preserve
bysort facility_id: egen min_rate = min(rate)
bysort facility_id: egen max_rate = max(rate)
collapse min_rate* max_rate*, by(facility_id)

sort min_rate
gen trim_1 = 1 if _n <= round(_N/100)
gen trim_5 = 1 if _n <= round(_N/20)
sort max_rate
replace trim_1 = 1 if _n > _N - round(_N/100)
replace trim_5 = 1 if _n > _N - round(_N/20)

drop min* max*
sort facility_id
save "$DATA/Temp/trim_facilities", replace
restore

merge m:1 facility_id using "$DATA/Temp/trim_facilities"
drop _merge

xtset facility_id

compress*
save "$DATA/ghgrp_di_intermediate", replace

// Preliminary regression macros
local dependent			rate
local independent  		p_spot
local controls		 	wells oil completions
local FE 				i.region_year
local conditions		if trim_5 != 1

// Recover average emission rate and price to vertically adjust graphs
qui sum `dependent' `conditions' 
local avg_rate `=r(mean)'
qui sum `independent' `conditions'
local avg_price `=r(mean)'



// 2 // Main Regressions

// 2.1 // Second-order FP
fp <`independent'>, dimension(2) powers(-2 -2) replace: xtreg `dependent' <`independent'> `controls' `FE' `conditions', fe vce(cluster facility)

mat B_2 = e(b) // Note: Graphing Requires extra steps to center because fixed effects shift it vertically - first recover center point of line (at average price), then adjust line to match average rate in sample
mat FP_2 =  e(fp_fp)
if FP_2[1,1] != FP_2[1,2] { // This calculates the fitted rate at sample average price (equation format depends on FP results - note, doesn't work for log)
	local rate_midpoint = B_2[1,`=colsof(B_2)'] + B_2[1,1]*(`avg_price'^`=FP_2[1,1]') + B_2[1,2]*(`avg_price'^`=FP_2[1,2]') 
}
if FP_2[1,1] == FP_2[1,2] {
	local rate_midpoint = B_2[1,`=colsof(B_2)'] + B_2[1,1]*(`avg_price'^`=FP_2[1,1]') + B_2[1,2]*log(`avg_price')*(`avg_price'^`=FP_2[1,1]')
}
local adjust `=`rate_midpoint' - `avg_rate'' // This is used to adjust the graph y-axis to center the regression line at the average emission rate
local original B_2[1,`=colsof(B_2)'] + B_2[1,1]*x^FP_2[1,1] + B_2[1,2]*log(x)*x^FP_2[1,2]

// Note:
// The previous ~10 lines are used to recover parameters used to adjust graphs for most robustness checks and add the dotted blue line, so do not comment them out (the rest of section 2 and each later section can be commented out to run robustness checks one group at a time)

/*
fp plot, residuals(none) plotopts(msize(small) mcolor(green)) lineopts(lwidth(medthick) lcolor(gs4)) ///
		ciopts(fcolor(gs13) lcolor(gs13)) ///
		graphregion(color(white)) ///
		xlabel(1(1)6, labsize(medlarge) format(%3.0fc)) xtitle("Price ($/Mcf)", size(large) margin(0 0 -2 2)) ///
		ytitle("Fitted Emission Rate", size(large) margin(-2 3 0 0)) ///
		ylabel(`adjust' "0" `=`adjust'+.005' ".005"  `=`adjust'+.01' ".01"  `=`adjust'+.015' ".015" `=`adjust'+.02' ".02", angle(vert) labsize(medlarge) glcolor(gs14) glwidth(vvthin) ) ///
		yscale(range(`=`adjust'-.003' `=`adjust'+.0201')) ///
		saving("Graphs/fp_main", replace) name(fp_main)
	
// 2.2 // Linear

fp <`independent'>, dimension(1) powers(1) replace: xtreg `dependent' <`independent'> `controls' `FE'  `conditions', fe vce(cluster facility)
mat B_0 = e(b)
local rate_midpoint_0 = B_0[1,`=colsof(B_0)'] + B_0[1,1]*(`avg_price')
local adjust_0 `=`rate_midpoint' - `rate_midpoint_0'' // This is used to adjust the linear regression line to match the same center point as the second-order FP
di `rate_midpoint_0'
di `adjust_0'


// 2.3 // First-order FP

fp <`independent'>, replace dimension(1) powers(): xtreg `dependent' <`independent'> `controls' `FE'  `conditions', fe vce(cluster facility)
mat B_1 = e(b)
mat FP_1 = e(fp_fp)
local rate_midpoint_1 = B_1[1,`=colsof(B_1)'] + B_1[1,1]*(`avg_price'^`=FP_1[1,1]')
local adjust_1 `=`rate_midpoint' - `rate_midpoint_1''
di `rate_midpoint_1'
di `adjust_1'


// 2.4 // Third-order FP

fp <`independent'>, replace dimension(3) powers(-.5 0): xtreg `dependent' <`independent'> `controls' `FE'  `conditions', fe vce(cluster facility)
mat B_3 = e(b)
mat FP_3 = e(fp_fp)
local rate_midpoint_3 = B_3[1,`=colsof(B_3)'] + B_3[1,1]*(`avg_price'^`=FP_3[1,1]') + B_3[1,2]*log(`avg_price') + B_3[1,3]*log(`avg_price')*log(`avg_price')
local adjust_3 `=`rate_midpoint' - `rate_midpoint_3''
di `rate_midpoint_3'
di `adjust_3'


// 2.5 // FP Comparison graph
qui sum p_spot
twoway  function y = `adjust_0' + B_0[1,`=colsof(B_0)'] + B_0[1,1]*x, lcolor(erose) lwidth(vthick) range(`r(min)' `r(max)') || ///
		function y = `adjust_1' + B_1[1,`=colsof(B_1)'] + B_1[1,1]*x^FP_1[1,1],  lcolor(edkblue) lwidth(thick) lpattern("-##") range(`r(min)' `r(max)') || ///
		function y = 			  B_2[1,`=colsof(B_2)'] + B_2[1,1]*x^FP_2[1,1] + B_2[1,2]*log(x)*x^FP_2[1,2], lcolor(gs3) lwidth(thick) range(`r(min)' `r(max)') || ///
		function y = `adjust_3' + B_3[1,`=colsof(B_3)'] + B_3[1,1]*x^FP_3[1,1] + B_3[1,2]*log(x) + B_3[1,3]*log(x)*log(x), lcolor(emidblue) lpattern("_##") lwidth(thick) range(`r(min)' `r(max)') ///	
		graphregion(color(white)) ///
		xlabel(1(1)6, labsize(medlarge) format(%3.0fc)) xtitle("Price ($/Mcf)", size(large) margin(0 0 -2 2)) ///
		ytitle(" ", size(large) margin(0 3 0 0)) ///
		yscale(range(`=`adjust'-.003' `=`adjust'+.0201')) ///
		ylabel(`adjust' "0" `=`adjust'+.005' ".005"  `=`adjust'+.01' ".01"  `=`adjust'+.015' ".015" `=`adjust'+.02' ".02", angle(vert) labsize(medlarge) glcolor(gs14) glwidth(vvthin) ) ///
		legend(order(1 "Linear" 2 "First-Order FP" 3 "Second-Order FP" 4 "Third-Order FP") size(medlarge) row(4) ring(0) pos(2) region( lwidth(thin) lcolor(gs14) margin( 1 1 1 1 ) )  ) ///
		saving("Graphs/fp_all", replace) name(fp_all)

graph combine "Graphs/fp_main.gph" "Graphs/fp_all.gph", imargin(2 2 2 2) iscale(1) ///
	graphregion(color(white)) ///
	xsize(10) ysize(4)

graph export "$DATA/Results/p8_fp_main.pdf", replace


// 2.6 // Regressions for table output

eststo clear

gen p_spot_half = p_spot^.5 // First-order FP

gen p_spot_n2 = p_spot^(-2) // Second-order FP
gen p_spot_n2_log = log(p_spot)*p_spot^(-2)

gen p_spot_nhalf = p_spot^(-.5) // Third-order FP
gen p_spot_log = log(p_spot)
gen p_spot_log_log = log(p_spot)*log(p_spot)

eststo: xtreg `dependent' p_spot `controls' i.region_year  `conditions', fe vce(cluster facility)
eststo: xtreg `dependent' p_spot_half `controls' i.region_year  `conditions', fe vce(cluster facility)
eststo: xtreg `dependent' p_spot_n2 p_spot_n2_log `controls' i.region_year  `conditions', fe vce(cluster facility)
eststo: xtreg `dependent' p_spot_nhalf p_spot_log p_spot_log_log `controls' i.region_year `conditions', fe vce(cluster facility)

esttab using "$DATA/Results/p8_main.tex", b(4) se(4) replace star(* 0.10 ** 0.05 *** 0.01) keep(p_spot* ) booktabs scalars(r2_w)
a
*/






// 3 // Robustness Checks 1 - Alternative standard errors
/*
eststo clear

// 3.1 // EHW
eststo: xtreg `dependent' `independent' `controls' i.region_year `conditions', fe vce(robust)

fp <`independent'>, dimension(2) powers(-2 -2) replace: xtreg `dependent' <`independent'> `controls' `FE'  `conditions', fe vce(robust)
mat B = e(b)
mat FP = e(fp_fp)
local rate_midpoint = B[1,`=colsof(B)'] + B[1,1]*(`avg_price'^`=FP[1,1]') + B[1,2]*log(`avg_price')*(`avg_price'^`=FP[1,1]')
local adjust `=`rate_midpoint' - `avg_rate'' // This is used to adjust the graph y-axis to center the regression line at the average emission rate
fp plot, residuals(none) plotopts(msize(small) mcolor(green)) lineopts(lwidth(medthick) lcolor(gs4)) ///
	ciopts(fcolor(gs13) lcolor(gs13)) ///
	graphregion(color(white)) ///
	xlabel(1(1)6, labsize(medlarge) format(%3.0fc)) xtitle("Price ($/Mcf)", size(large) margin(0 0 -2 2)) ///
	ytitle("Fitted Emission Rate", size(large) margin(-2 3 0 0)) ///
	yscale(range(`=`adjust'-.003' `=`adjust'+.0231')) ///
	ylabel(`adjust' "0" `=`adjust'+.005' ".005"  `=`adjust'+.01' ".01"  `=`adjust'+.015' ".015" `=`adjust'+.02' ".02", angle(vert) labsize(medlarge) glcolor(gs14) glwidth(vvthin) ) ///
	title("Heteroskedasticity-Consistent (EHW)", size(large)) ///
	name(fp_se_hwe) ///
	saving(Graphs/fp_se_hwe, replace)



// 3.2 // Clustering at parent firm level
eststo: xtreg `dependent' `independent' `controls' i.region_year `conditions', fe vce(cluster facility)

fp <`independent'>, dimension(2) powers(-2 -2) replace: xtreg `dependent' <`independent'> `controls' `FE' `conditions', fe vce(cluster parent)
mat B = e(b)
mat FP = e(fp_fp)
local rate_midpoint = B[1,`=colsof(B)'] + B[1,1]*(`avg_price'^`=FP[1,1]') + B[1,2]*log(`avg_price')*(`avg_price'^`=FP[1,1]')
local adjust `=`rate_midpoint' - `avg_rate'' // This is used to adjust the graph y-axis to center the regression line at the average emission rate
fp plot, residuals(none) plotopts(msize(small) mcolor(green)) lineopts(lwidth(medthick) lcolor(gs4)) ///
	ciopts(fcolor(gs13) lcolor(gs13)) ///
	graphregion(color(white)) ///
	xlabel(1(1)6, labsize(medlarge) format(%3.0fc)) xtitle("Price ($/Mcf)", size(large) margin(0 0 -2 2)) ///
	ytitle("Fitted Emission Rate", size(large) margin(-2 3 0 0)) ///
	yscale(range(`=`adjust'-.003' `=`adjust'+.0231')) ///
	ylabel(`adjust' "0" `=`adjust'+.005' ".005"  `=`adjust'+.01' ".01"  `=`adjust'+.015' ".015" `=`adjust'+.02' ".02", angle(vert) labsize(medlarge) glcolor(gs14) glwidth(vvthin) ) ///
	title("SE Clustered at Parent Firm Level", size(large)) ///
	name(fp_se_parent) ///
	saving(Graphs/fp_se_parent, replace)

// 3.3 // Clustering at basin level
eststo: xtreg `dependent' `independent' `controls' i.region_year `conditions', fe vce(cluster basin)
fp <`independent'>, dimension(2) powers(-2 -2) replace: xtreg `dependent' <`independent'> `controls' `FE' `conditions', fe vce(cluster basin)
mat B = e(b)
mat FP = e(fp_fp)
local rate_midpoint = B[1,`=colsof(B)'] + B[1,1]*(`avg_price'^`=FP[1,1]') + B[1,2]*log(`avg_price')*(`avg_price'^`=FP[1,1]')
local adjust `=`rate_midpoint' - `avg_rate'' // This is used to adjust the graph y-axis to center the regression line at the average emission rate
fp plot, residuals(none) plotopts(msize(small) mcolor(green)) lineopts(lwidth(medthick) lcolor(gs4)) ///
	ciopts(fcolor(gs13) lcolor(gs13)) ///
	graphregion(color(white)) ///
	xlabel(1(1)6, labsize(medlarge) format(%3.0fc)) xtitle("Price ($/Mcf)", size(large) margin(0 0 -2 2)) ///
	ytitle("Fitted Emission Rate", size(large) margin(-2 3 0 0)) ///
	yscale(range(`=`adjust'-.003' `=`adjust'+.0231')) ///
	ylabel(`adjust' "0" `=`adjust'+.005' ".005"  `=`adjust'+.01' ".01"  `=`adjust'+.015' ".015" `=`adjust'+.02' ".02", angle(vert) labsize(medlarge) glcolor(gs14) glwidth(vvthin) ) ///
	title("SE Clustered at Basin Level", size(large)) ///
	saving(Graphs/fp_se_basin, replace) name(fp_se_basin)

// 3.4 // Clustering at year level
eststo: xtreg `dependent' `independent' `controls' i.region_year `conditions', nonest fe vce(cluster year)

fp <`independent'>, dimension(2) powers(-2 -2) replace: xtreg `dependent' <`independent'> `controls' `FE'  `conditions', nonest fe vce(cluster year)
mat B = e(b)
mat FP = e(fp_fp)
local rate_midpoint = B[1,`=colsof(B)'] + B[1,1]*(`avg_price'^`=FP[1,1]') + B[1,2]*log(`avg_price')*(`avg_price'^`=FP[1,1]')
local adjust `=`rate_midpoint' - `avg_rate'' // This is used to adjust the graph y-axis to center the regression line at the average emission rate
fp plot, residuals(none) plotopts(msize(small) mcolor(green)) lineopts(lwidth(medthick) lcolor(gs4)) ///
	ciopts(fcolor(gs13) lcolor(gs13)) ///
	graphregion(color(white)) ///
	xlabel(1(1)6, labsize(medlarge) format(%3.0fc)) xtitle("Price ($/Mcf)", size(large) margin(0 0 -2 2)) ///
	ytitle("Fitted Emission Rate", size(large) margin(-2 3 0 0)) ///
	yscale(range(`=`adjust'-.003' `=`adjust'+.0231')) ///
	ylabel(`adjust' "0" `=`adjust'+.005' ".005"  `=`adjust'+.01' ".01"  `=`adjust'+.015' ".015" `=`adjust'+.02' ".02", angle(vert) labsize(medlarge) glcolor(gs14) glwidth(vvthin) ) ///
	title("SE Clustered at Year Level", size(large)) ///
	saving(Graphs/fp_se_year, replace) name(fp_se_year)
	
// 3.5 // Trimming emission rates at 1%

use "$DATA/ghgrp_di_intermediate", clear
local conditions if trim_1 != 1
eststo: xtreg `dependent' `independent' `controls' i.region_year `conditions', fe vce(cluster facility)

fp <`independent'>, dimension(2) powers() replace: xtreg `dependent' <`independent'> `controls' `FE'  `conditions', fe vce(cluster facility)
mat B = e(b)
mat FP = e(fp_fp)
local rate_midpoint_R = B[1,`=colsof(B)'] + B[1,1]*(`avg_price'^FP[1,1]) + B[1,2]*(`avg_price'^FP[1,2])*log(`avg_price') // Note - this one has repeated power, so second is multiplied by log
local adjust_ax `=`rate_midpoint_R' - `avg_rate''
local adjust_orig `=`rate_midpoint_R' - `rate_midpoint''

fp plot, residuals(none) plotopts(msize(small) mcolor(green)) lineopts(lwidth(medthick) lcolor(gs4)) ///
	ciopts(fcolor(gs13) lcolor(gs13) ) ///
	graphregion(color(white)) name(fp_trim_1) ///
	xlabel(1(1)6, labsize(medlarge) format(%3.0fc)) xtitle("Price ($/Mcf)", size(large) margin(0 0 -2 2)) ///
	ytitle("Fitted Emission Rate", size(large) margin(0 2 0 0)) ///
	yscale(range(`=`adjust'-.003' `=`adjust'+.0231')) ///
	ylabel(`adjust' "0" `=`adjust'+.005' ".005"  `=`adjust'+.01' ".01"  `=`adjust'+.015' ".015" `=`adjust'+.02' ".02", angle(vert) labsize(medlarge) glcolor(gs14) glwidth(vvthin) ) ///
	title("Trimming Emission Rates at 1%", size(large))

qui sum p_spot
graph addplot function y = `adjust_orig' + `original', lcolor(ebblue) lpattern("_##") lwidth(thick) range(`r(min)' `r(max)') ///
	ylabel(`=`adjust_ax'-.04' "-.04"  `=`adjust_ax'-.02' "-.02" `adjust_ax' "0" `=`adjust_ax'+.02' ".02" `=`adjust_ax'+.04' ".04"  `=`adjust_ax'+.06' ".06" , angle(vert) labsize(medlarge) nogrid glcolor(white) glwidth(vvthin) ) ///
	ymlabel(`=`adjust_ax'-.04'(.005)`=`adjust_ax'+.06', nolabel tlcolor(white) grid glcolor(gs14) glwidth(vvthin)) ///
	yscale(range(`=`adjust_ax'-0.0001' `=`adjust_ax'+0.0601')) ///
	saving(Graphs/fp_trim_1, replace)

local conditions if trim_5 != 1 // Reset trim

// 3.6 // Combine and export results
graph combine Graphs/fp_se_hwe.gph Graphs/fp_se_parent.gph Graphs/fp_se_basin.gph Graphs/fp_se_year.gph Graphs/fp_trim_1.gph, ///
	col(2) imargin(2 4 6 2) iscale(.4) ///
	graphregion(color(white)) ///
	xsize(8.5) ysize(10)

graph export "$DATA/Results/p8_robust_1.pdf", replace

esttab using "$DATA/Results/p8_robust_1.tex", b(4) se(4) replace star(* 0.10 ** 0.05 *** 0.01) keep(p_spot) booktabs scalars(r2_w)
a
*/








// 4 // Robustness Checks 2 - Misc identification questions
/*
eststo clear

// 4.1 // State Regs

use "$DATA/ghgrp_di_intermediate", clear

drop if state_regs == 1

preserve
bysort facility_id: egen min_rate = min(rate)
bysort facility_id: egen max_rate = max(rate)
collapse min_rate max_rate, by(facility_id)
sort min_rate
gen trim_5 = 1 if _n <= round(_N/20)
sort max_rate
replace trim_5 = 1 if _n > _N - round(_N/20)
drop min max
sort facility_id
save "$DATA/Temp/trim_facilities", replace
restore
merge m:1 facility_id using "$DATA/Temp/trim_facilities"
drop _merge

eststo: xtreg `dependent' `independent' `controls' i.region_year `conditions', fe vce(cluster facility)

fp <`independent'>, dimension(2) powers(): xtreg `dependent' <`independent'> `controls' `FE' `conditions', fe vce(cluster facility)
mat B = e(b)
mat FP = e(fp_fp)
local rate_midpoint_R = B[1,`=colsof(B)'] + B[1,1]*(`avg_price'^FP[1,1]) + B[1,2]*(`avg_price'^FP[1,2])*log(`avg_price')
local adjust_ax `=`rate_midpoint_R' - `avg_rate''
local adjust_orig `=`rate_midpoint_R' - `rate_midpoint''

fp plot, residuals(none) plotopts(msize(small) mcolor(green)) lineopts(lwidth(medthick) lcolor(gs4)) ///
	ciopts(fcolor(gs13) lcolor(gs13) ) ///
	graphregion(color(white)) name(fp_regs) ///
	xlabel(1(1)6, labsize(medlarge) format(%3.0fc)) xtitle("Price ($/Mcf)", size(large) margin(0 0 -2 2)) ///
	ytitle("Fitted Emission Rate", size(large) margin(0 2 0 0)) ///
	title("Trimming Based on State Regulations", size(large)) 
	
qui sum p_spot
graph addplot function y = `adjust_orig' + `original', lcolor(ebblue) lpattern("_##") lwidth(thick) range(`r(min)' `r(max)') ///
	ylabel(`=`adjust_ax'-.005' "-.005" `adjust_ax' "0" `=`adjust_ax'+.005' ".005"  `=`adjust_ax'+.01' ".01"  `=`adjust_ax'+.015' ".015" `=`adjust_ax'+.02' ".02", angle(vert) labsize(medlarge) glcolor(gs14) glwidth(vvthin) ) ///
	yscale(range(`=`adjust_ax'-0.0071' `=`adjust_ax'+0.022')) ///
	saving(Graphs/fp_regs, replace)

// 4.2 // Controlling for completions in the previous year

use "$DATA/ghgrp_di_intermediate", clear

drop if missing(ghgrp_comp_lastyr)
preserve
bysort facility_id: egen min_rate = min(rate)
bysort facility_id: egen max_rate = max(rate)
collapse min_rate max_rate, by(facility_id)
sort min_rate
gen trim_5 = 1 if _n <= round(_N/20)
sort max_rate
replace trim_5 = 1 if _n > _N - round(_N/20)
drop min max
sort facility_id
save "$DATA/Temp/trim_facilities", replace
restore
merge m:1 facility_id using "$DATA/Temp/trim_facilities"
drop _merge

local controls wells oil completions ghgrp_comp_lastyr
eststo: xtreg `dependent' `independent' `controls' i.region_year `conditions', fe vce(cluster facility)

fp <`independent'>, dimension(2) powers(): xtreg `dependent' <`independent'> `controls' `FE'  `conditions', fe vce(cluster facility)
mat B = e(b)
mat FP = e(fp_fp)
local rate_midpoint_R = B[1,`=colsof(B)'] + B[1,1]*(`avg_price'^FP[1,1]) + B[1,2]*(`avg_price'^FP[1,2]) 
local adjust_ax `=`rate_midpoint_R' - `avg_rate''
local adjust_orig `=`rate_midpoint_R' - `rate_midpoint''

fp plot, residuals(none) plotopts(msize(small) mcolor(green)) lineopts(lwidth(medthick) lcolor(gs4)) ///
	ciopts(fcolor(gs13) lcolor(gs13) ) ///
	graphregion(color(white)) name(fp_comp_lastyr) ///
	xlabel(1(1)6, labsize(medlarge) format(%3.0fc)) xtitle("Price ($/Mcf)", size(large) margin(0 0 -2 2)) ///
	ytitle("Fitted Emission Rate", size(large) margin(0 2 0 0)) ///
	title("Controlling for Completions in Previous Year", size(large)) 
	
qui sum p_spot
graph addplot function y = `adjust_orig' + `original', lcolor(ebblue) lpattern("_##") lwidth(thick) range(`r(min)' `r(max)') ///
	ylabel( `=`adjust_ax'-.005' "-.005" `adjust_ax' "0" `=`adjust_ax'+.005' ".005"  `=`adjust_ax'+.01' ".01"  `=`adjust_ax'+.015' ".015" `=`adjust_ax'+.02' ".02" , angle(vert) labsize(medlarge) glcolor(gs14) glwidth(vvthin) ) ///
	yscale(range(`=`adjust_ax'-0.0071' `=`adjust_ax'+0.022')) ///
	saving(Graphs/fp_comp_lastyr, replace)

local controls wells oil completions

// 4.3 // Year FE 

use "$DATA/ghgrp_di_intermediate", clear
local FE i.year

eststo: xtreg `dependent' `independent' `controls' i.year `conditions', fe vce(cluster facility)

fp <`independent'>, dimension(2) powers() replace: xtreg `dependent' <`independent'> `controls' `FE'  `conditions', fe vce(cluster facility)
mat B = e(b)
mat FP = e(fp_fp)
local rate_midpoint_R = B[1,`=colsof(B)'] + B[1,1]*(`avg_price'^FP[1,1]) + B[1,2]*(`avg_price'^FP[1,2])*log(`avg_price') // Note - this one has repeated power, so second is multiplied by log
local adjust_ax `=`rate_midpoint_R' - `avg_rate''
local adjust_orig `=`rate_midpoint_R' - `rate_midpoint''

fp plot, residuals(none) plotopts(msize(small) mcolor(green)) lineopts(lwidth(medthick) lcolor(gs4)) ///
	ciopts(fcolor(gs13) lcolor(gs13) ) ///
	graphregion(color(white)) name(fp_year) ///
	xlabel(1(1)6, labsize(medlarge) format(%3.0fc)) xtitle("Price ($/Mcf)", size(large) margin(0 0 -2 2)) ///
	ytitle("Fitted Emission Rate", size(large) margin(0 2 0 0)) ///
	title("Year FE", size(large))

qui sum p_spot
graph addplot function y = `adjust_orig' + `original', lcolor(ebblue) lpattern("_##") lwidth(thick) range(`r(min)' `r(max)') ///
	ylabel( `=`adjust_ax'-.005' "-.005" `adjust_ax' "0" `=`adjust_ax'+.005' ".005"  `=`adjust_ax'+.01' ".01"  `=`adjust_ax'+.015' ".015" `=`adjust_ax'+.02' ".02" , angle(vert) labsize(medlarge) glcolor(gs14) glwidth(vvthin) ) ///
	yscale(range(`=`adjust_ax'-0.0071' `=`adjust_ax'+0.022')) ///
	saving(Graphs/fp_year_fe, replace)

local FE i.region_year i.facility_id


// 4.4 // Basin-Year FE

use "$DATA/ghgrp_di_intermediate", clear
local FE i.basin_year

eststo: xtreg `dependent' `independent' `controls' i.basin_year `conditions', fe vce(cluster facility)

fp <`independent'>, dimension(2) powers() replace: xtreg `dependent' <`independent'> `controls' `FE'  `conditions', fe vce(cluster facility)
mat B = e(b)
mat FP = e(fp_fp)
local rate_midpoint_R = B[1,`=colsof(B)'] + B[1,1]*(`avg_price'^FP[1,1]) + B[1,2]*(`avg_price'^FP[1,2])*log(`avg_price') // Note - this one has repeated power, so second is multiplied by log
local adjust_ax `=`rate_midpoint_R' - `avg_rate''
local adjust_orig `=`rate_midpoint_R' - `rate_midpoint''

fp plot, residuals(none) plotopts(msize(small) mcolor(green)) lineopts(lwidth(medthick) lcolor(gs4)) ///
	ciopts(fcolor(gs13) lcolor(gs13) ) ///
	graphregion(color(white)) name(fp_basin_year) ///
	xlabel(1(1)6, labsize(medlarge) format(%3.0fc)) xtitle("Price ($/Mcf)", size(large) margin(0 0 -2 2)) ///
	ytitle("Fitted Emission Rate", size(large) margin(0 2 0 0)) ///
	title("Basin-Year FE", size(large)) ///
	
qui sum p_spot
graph addplot function y = `adjust_orig' + `original', lcolor(ebblue) lpattern("_##") lwidth(thick) range(`r(min)' `r(max)') ///
	ylabel( `=`adjust_ax'-.005' "-.005" `adjust_ax' "0" `=`adjust_ax'+.005' ".005"  `=`adjust_ax'+.01' ".01"  `=`adjust_ax'+.015' ".015" `=`adjust_ax'+.02' ".02" , angle(vert) labsize(medlarge) glcolor(gs14) glwidth(vvthin) ) ///
	yscale(range(`=`adjust_ax'-0.0071' `=`adjust_ax'+0.022')) ///
	saving(Graphs/fp_basin_year_fe, replace)

local FE i.region_year i.facility_id


// 4.5 // Excluding mountain region (more code because observations must be added or dropped before trimming at 5%)

use "$DATA/ghgrp_di_intermediate", clear

drop if region == 4

preserve
bysort facility_id: egen min_rate = min(rate)
bysort facility_id: egen max_rate = max(rate)
collapse min_rate max_rate, by(facility_id)
sort min_rate
gen trim_5 = 1 if _n <= round(_N/20)
sort max_rate
replace trim_5 = 1 if _n > _N - round(_N/20)
drop min max
sort facility_id
save "$DATA/Temp/trim_facilities", replace
restore
merge m:1 facility_id using "$DATA/Temp/trim_facilities"
drop _merge

eststo: xtreg `dependent' `independent' `controls' i.region_year `conditions', fe vce(cluster facility)

fp <`independent'>, dimension(2) powers(): xtreg `dependent' <`independent'> `controls' `FE' `conditions', fe vce(cluster facility)
mat B = e(b)
mat FP = e(fp_fp)
local rate_midpoint_R = B[1,`=colsof(B)'] + B[1,1]*(`avg_price'^FP[1,1]) + B[1,2]*(`avg_price'^FP[1,2])*log(`avg_price') // Note - this one has repeated power, so second is multiplied by log
local adjust_ax `=`rate_midpoint_R' - `avg_rate''
local adjust_orig `=`rate_midpoint_R' - `rate_midpoint''

fp plot, residuals(none) plotopts(msize(small) mcolor(green)) lineopts(lwidth(medthick) lcolor(gs4)) ///
	ciopts(fcolor(gs13) lcolor(gs13) ) ///
	graphregion(color(white)) name(fp_no_mountain) ///
	xlabel(1(1)6, labsize(medlarge) format(%3.0fc)) xtitle("Price ($/Mcf)", size(large) margin(0 0 -2 2)) ///
	ytitle("Fitted Emission Rate", size(large) margin(0 2 0 0)) ///
	title("Excluding Mountain Region", size(large))

qui sum p_spot
graph addplot function y = `adjust_orig' + `original', lcolor(ebblue) lpattern("_##") lwidth(thick) range(`r(min)' `r(max)') ///
	ylabel( `=`adjust_ax'-.005' "-.005" `adjust_ax' "0" `=`adjust_ax'+.005' ".005"  `=`adjust_ax'+.01' ".01"  `=`adjust_ax'+.015' ".015" `=`adjust_ax'+.02' ".02" , angle(vert) labsize(medlarge) glcolor(gs14) glwidth(vvthin) ) ///
	yscale(range(`=`adjust_ax'-0.0071' `=`adjust_ax'+0.022')) ///
	saving(Graphs/fp_no_mountain, replace)


// 4.6 // Dropping 2016

use "$DATA/ghgrp_di_intermediate", clear

drop if year == 2016
preserve
bysort facility_id: egen min_rate = min(rate)
bysort facility_id: egen max_rate = max(rate)
collapse min_rate max_rate, by(facility_id)
sort min_rate
gen trim_5 = 1 if _n <= round(_N/20)
sort max_rate
replace trim_5 = 1 if _n > _N - round(_N/20)
drop min max
sort facility_id
save "$DATA/Temp/trim_facilities", replace
restore
merge m:1 facility_id using "$DATA/Temp/trim_facilities"
drop _merge

eststo: xtreg `dependent' `independent' `controls' i.region_year `conditions', fe vce(cluster facility)

fp <`independent'>, dimension(2) powers(): xtreg `dependent' <`independent'> `controls' `FE'  `conditions', fe vce(cluster facility)
mat B = e(b)
mat FP = e(fp_fp)
local rate_midpoint_R = B[1,`=colsof(B)'] + B[1,1]*(`avg_price'^FP[1,1]) + B[1,2]*(`avg_price'^FP[1,2])*log(`avg_price') // Note - this one has repeated power, so second is multiplied by log
local adjust_ax `=`rate_midpoint_R' - `avg_rate''
local adjust_orig `=`rate_midpoint_R' - `rate_midpoint''

fp plot, residuals(none) plotopts(msize(small) mcolor(green)) lineopts(lwidth(medthick) lcolor(gs4)) ///
	ciopts(fcolor(gs13) lcolor(gs13) ) ///
	graphregion(color(white)) name(fp_no_2016) ///
	xlabel(1(1)6, labsize(medlarge) format(%3.0fc)) xtitle("Price ($/Mcf)", size(large) margin(0 0 -2 2)) ///
	ytitle("Fitted Emission Rate", size(large) margin(0 2 0 0)) ///
	title("Excluding 2016", size(large)) 

qui sum p_spot
graph addplot function y = `adjust_orig' + `original', lcolor(ebblue) lpattern("_##") lwidth(thick) range(`r(min)' `r(max)') ///
	ylabel( `=`adjust_ax'-.005' "-.005" `adjust_ax' "0" `=`adjust_ax'+.005' ".005"  `=`adjust_ax'+.01' ".01"  `=`adjust_ax'+.015' ".015" `=`adjust_ax'+.02' ".02" , angle(vert) labsize(medlarge) glcolor(gs14) glwidth(vvthin) ) ///
	yscale(range(`=`adjust_ax'-0.0071' `=`adjust_ax'+0.022')) ///
	saving(Graphs/fp_no_2016, replace)

// 4.7 // Combining and exporting
graph combine Graphs/fp_regs.gph Graphs/fp_comp_lastyr.gph Graphs/fp_year_fe.gph Graphs/fp_basin_year_fe.gph Graphs/fp_no_mountain.gph Graphs/fp_no_2016.gph, ///
	col(2) imargin(2 4 6 2) iscale(.4) ///
	graphregion(color(white)) ///
	xsize(8.5) ysize(10)
graph export "$DATA/Results/p8_robust_2.pdf", replace

esttab using "Results/p8_robust_2.tex", b(4) se(4) replace star(* 0.10 ** 0.05 *** 0.01) keep(p_spot) booktabs scalars(r2_w)
a
*/






// 5 // Negative Binomial Model
/*
eststo clear

rename ch4_normal ch4
local FE i.region_year i.facility_id

eststo: fp <`independent'>, dimension(1) powers(1) replace: nbreg ch4 <`independent'> `controls' `FE' `conditions', exposure(gas) vce(cluster facility) // Linear
mat NB_0 = e(b)
eststo: fp <`independent'>, dimension(1) powers() replace: nbreg ch4 <`independent'> `controls' `FE' `conditions', exposure(gas) vce(cluster facility) // 1st-degree FP
mat NB_1 = e(b)
eststo: fp <`independent'>, dimension(2) powers() replace: nbreg ch4 <`independent'> `controls' `FE' `conditions', exposure(gas) vce(cluster facility) // 2nd-degree FP
mat NB_2 = e(b)
eststo: fp <`independent'>, dimension(3) powers() replace: nbreg ch4 <`independent'> `controls' `FE' `conditions', exposure(gas) vce(cluster facility) // 3rd-degree FP
mat NB_3 = e(b)

esttab using "Results\p8_nb.tex", b(4) se(4) replace star(* 0.10 ** 0.05 *** 0.01) keep(p_spot*) booktabs ar2 transform(p_spot (exp(@)-1) exp(@))


gen log_rate = log(rate)
sum log_rate if trim_5 != 1 // -5.5
// Approximate the graph adjustment in this one to center all on -5.5. Note that this only changes intercept, not slope, which is the important part
	
fp <`independent'>, dimension(2) powers(-2 -1) replace: nbreg ch4 <`independent'> `controls' `FE' `conditions', exposure(gas) vce(cluster facility) // For graph
fp plot, residuals(none) plotopts(msize(small) mcolor(green)) lineopts(lwidth(medthick) lcolor(gs4)) ///
	ciopts(fcolor(gs13) lcolor(gs13)) ///
	graphregion(color(white)) ///
	xlabel(1(1)6, labsize(medlarge) format(%3.0fc)) xtitle("Price ($/Mcf)", size(large) margin(0 0 -2 2)) ///
	ytitle("log(Fitted Emission Rate)", size(large) margin(-2 3 0 0)) ///
	ylabel(-3.5 "-2" -5.5 "-4" -7.5 "-6" -9.5 "-8", labsize(medlarge) glcolor(gs14) glwidth(vvthin) ) ///
	yscale(range(-3.5 -9.5))  ///
	saving(Graphs/fp_nb, replace)

qui sum p_spot
twoway  function y = NB_0[1,`=colsof(NB_0)-1'] + NB_0[1,1]*x, lcolor(erose) lwidth(vthick) range(`r(min)' `r(max)') || ///
		function y = NB_1[1,`=colsof(NB_1)-1'] + NB_1[1,1]*x^(.5),  lcolor(edkblue) lwidth(thick) lpattern(dash) range(`r(min)' `r(max)') || ///
		function y = NB_2[1,`=colsof(NB_2)-1'] + NB_2[1,1]*x^(-2) + NB_2[1,2]*x^(-1), lcolor(gs3) lwidth(thick) range(`r(min)' `r(max)') || ///
		function y = NB_3[1,`=colsof(NB_3)-1'] + NB_3[1,1]*x^(-2) + NB_3[1,2]*x^(-2)*log(x) + NB_3[1,3]*x^(3), lcolor(emidblue) lpattern(longdash) lwidth(thick) range(`r(min)' `r(max)') ///	
		graphregion(color(white)) ///
		xlabel(1(1)6, labsize(medlarge) format(%3.0fc)) xtitle("Price ($/Mcf)", size(large) margin(0 0 -2 2)) ///
		ytitle(" ", size(large) margin(0 3 0 0)) ///
		ylabel(-3.5 "-2" -5.5 "-4" -7.5 "-6" -9.5 "-8", labsize(medlarge) glcolor(gs14) glwidth(vvthin) ) ///
		yscale(range(-3.5 -9.5))  ///
		legend(order(1 "Linear" 2 "First-Order FP" 3 "Second-Order FP" 4 "Third-Order FP") size(medlarge) row(4) ring(0) pos(2) region( lwidth(thin) lcolor(gs14) margin( 1 1 1 1 ) )  ) ///
		saving(Graphs/fp_nb_all, replace)

graph combine Graphs/fp_nb.gph Graphs/fp_nb_all.gph, imargin(2 2 2 2) iscale(1) ///
	graphregion(color(white)) ///
	xsize(10) ysize(4)

graph export "$DATA/Results/p8_nb.pdf", replace

local FE i.region_year
a
*/





// 6 // Abatement mechanisms
/*
eststo clear

use "$DATA/ghgrp_di", clear

// Merging prices (proprietary data - see ReadMe for instructions on how to incorporate)
merge m:1 facility_id year using "$DATA\Prices\ghgrp_prices"
drop if _m == 2
drop _m

xtset facility_id

// Dropping facilities with operations in either Colorado or ND
gen trim_co_nd = 1 if colorado != 0 | north_dakota != 0
drop if trim_co_nd == 1

eststo clear
foreach var of varlist pneum_pumps pd_lo pd_inter pd_hi {
	preserve
	drop if missing(`var')
	bysort facility_id: egen min_rate = min(`var')
	bysort facility_id: egen max_rate = max(`var')
	collapse min_rate* max_rate* gas, by(facility_id)
	gsort min_rate -gas // Make sure to drop the same observations that have 0 each time - choose least plausible facilities
	gen trim_5 = 1 if _n <= round(_N/20)
	sort max_rate
	replace trim_5 = 1 if _n > _N - round(_N/20)
	drop min* max*
	sort facility_id
	save "$DATA/Temp/trim_facilities", replace
	restore
	merge m:1 facility_id using "$DATA/Temp/trim_facilities"
	drop _merge	
	eststo: xtreg `var' p_spot wells oil completions i.region_year `conditions', fe vce(cluster facility)
	drop trim_5
}
esttab, se replace star(* 0.10 ** 0.05 *** 0.01) keep(p_spot) scalars(r2_w)
esttab using "$DATA/Results/p8_mech.tex", se replace star(* 0.10 ** 0.05 *** 0.01) keep(p_spot) booktabs scalars(r2_w)
a
*/



// 7 // Other Sources
/*
eststo clear

// 7.1 // Completions
use "$DATA/ghgrp_di", clear

// Merging prices (proprietary data - see ReadMe for instructions on how to incorporate)
merge m:1 facility_id year using "$DATA/Prices/ghgrp_prices"
drop if _m == 2
drop _m

xtset facility_id

// Dropping facilities with operations in either Colorado or ND
gen trim_co_nd = 1 if colorado != 0 | north_dakota != 0
drop if trim_co_nd == 1

// Trimming emission rate outliers at the facility level
preserve
drop if missing(rate_comp)
bysort facility_id: egen min_rate = min(rate_comp)
bysort facility_id: egen max_rate = max(rate_comp)
collapse min_rate* max_rate* comp_work, by(facility_id)

sort min_rate comp_work // Need to ensure the same facilities are trimmed each time
gen trim_1 = 1 if _n <= round(_N/100)
gen trim_5 = 1 if _n <= round(_N/20)
sort max_rate
replace trim_1 = 1 if _n > _N - round(_N/100)
replace trim_5 = 1 if _n > _N - round(_N/20)

drop min* max*
sort facility_id
save "$DATA/Temp/trim_facilities", replace
restore

merge m:1 facility_id using "$DATA/Temp/trim_facilities"
drop _merge

// Preliminary regression macros
local dependent			rate_comp
local independent  		p_spot
local controls		 	wells oil
local FE 				i.region_year
local conditions		if trim_5 != 1

sum rate_comp if trim_5 != 1 // Mean 324 - Manually adjust y axis to center regression line here

eststo: xtreg `dependent' `independent' `controls' i.region_year `conditions', fe vce(cluster facility)

eststo: fp <`independent'>, dimension(2) powers() replace: xtreg `dependent' <`independent'> `controls' `FE' `conditions' , fe vce(cluster facility) 

fp plot, residuals(none) plotopts(msize(small) mcolor(green)) lineopts(lwidth(medthick) lcolor(gs4)) ///
	ciopts(fcolor(gs13) lcolor(gs13)) ///
	graphregion(color(white)) ///
	xlabel(1(1)6, labsize(medlarge) format(%6.0fc)) xtitle("Price ($/Mcf)", size(large) margin(0 0 -2 2)) ///
	ytitle("Predicted Emission Rate", size(large) margin(0 2 0 0)) ///
	title("Completions", size(large)) ///
	yscale(range()) ylabel(-1450 "-1000" -450 "0" 550 "1000" 1550 "2000", labsize(medlarge) glcolor(gs14) glwidth(vvthin) ) ///
	saving(Graphs/fp_comp, replace)


// 7.2 // Associated Gas

use "$DATA/ghgrp_di", clear

// Merging prices (proprietary data - see ReadMe for instructions on how to incorporate)
merge m:1 facility_id year using "$DATA/Prices/ghgrp_prices"
drop if _m == 2
drop _m

xtset facility_id

// Dropping facilities with operations in either Colorado or ND
gen trim_co_nd = 1 if colorado != 0 | north_dakota != 0
drop if trim_co_nd == 1

// Dropping facilities whose DI production values differ substantially from their GHGRP production values by more than 20%
gen prod_ratio = prod_mcf/gas if year >= 2015
bysort facility_id: egen mean_prod_ratio = mean(prod_ratio)
gen trim_prod = 1 if mean_prod_ratio < .8 | mean_prod_ratio > 1.25
drop prod_ratio mean_prod_ratio
drop if trim_prod == 1

// Trimming emission rate outliers at the facility level
preserve
drop if missing(rate_assoc)
bysort facility_id: egen min_rate = min(rate_assoc)
bysort facility_id: egen max_rate = max(rate_assoc)
collapse min_rate* max_rate*, by(facility_id)
sort min_rate
gen trim_1 = 1 if _n <= round(_N/100)
gen trim_5 = 1 if _n <= round(_N/20)
sort max_rate
replace trim_1 = 1 if _n > _N - round(_N/100)
replace trim_5 = 1 if _n > _N - round(_N/20)
drop min* max*
sort facility_id
save "$DATA/Temp/trim_facilities", replace
restore

merge m:1 facility_id using "$DATA/Temp/trim_facilities"
drop _merge

// Preliminary regression macros
local dependent			rate_assoc
local independent  		p_spot
local controls		 	wells oil completions
local FE 				i.region_year 
local conditions		if trim_5 != 1

eststo: xtreg `dependent' `independent' `controls' i.region_year `conditions', fe vce(cluster facility)

sum rate_assoc p_spot if trim_5 != 1 // Mean .014 - Manually adjust y axis to center regression line here

eststo: fp <`independent'>, dimension(2) powers() replace: xtreg `dependent' <`independent'> `controls' `FE' `conditions' , fe vce(cluster facility) 

fp plot, residuals(none) plotopts(msize(small) mcolor(green)) lineopts(lwidth(medthick) lcolor(gs4)) ///
	ciopts(fcolor(gs13) lcolor(gs13)) ///
	graphregion(color(white)) ///
	xlabel(1(1)6, labsize(medlarge) format(%6.0fc)) xtitle("Price ($/Mcf)", size(large) margin(0 0 -2 2)) ///
	ytitle("Predicted Emission Rate", size(large) margin(0 2 0 0)) ///
	title("Associated Gas", size(large)) ///
	yscale(range()) ylabel(-.06 "-.05" -.01 "0" .04 ".05" .09 ".1", labsize(medlarge) glcolor(gs14) glwidth(vvthin) ) ///
	saving(Graphs/fp_assoc, replace)
	
graph combine Graphs/fp_comp.gph Graphs/fp_assoc.gph, imargin(2 4 2 2) iscale(1) ///
	graphregion(color(white)) ///
	xsize(10) ysize(4)
graph export "$DATA/Results/p8_other_sources.pdf", replace

esttab using "Results/p8_sources.tex", se replace star(* 0.10 ** 0.05 *** 0.01) keep(p_spot*) booktabs scalars(r2_w)
a
*/



// 8 // Controlling for facility size
/*
eststo clear

// 8.1 // Weighted Regression

use "$DATA/ghgrp_di_intermediate", clear
drop if trim_5 == 1
eststo: xtreg `dependent' `independent' `controls' i.region_year [aweight=mean_gas] `conditions', fe vce(cluster facility)

eststo: fp <`independent'>, dimension(2) powers(-2 -.5) replace: xtreg `dependent' <`independent'> `controls' `FE' [aweight=mean_gas] `conditions', fe vce(cluster facility)
mat B = e(b)
mat FP = e(fp_fp)
local rate_midpoint_R = B[1,`=colsof(B)'] + B[1,1]*(`avg_price'^FP[1,1]) + B[1,2]*(`avg_price'^FP[1,2])
local adjust_ax `=`rate_midpoint_R' - `avg_rate''
local adjust_orig `=`rate_midpoint_R' - `rate_midpoint''

fp plot, residuals(none) plotopts(msize(small) mcolor(green)) lineopts(lwidth(medthick) lcolor(gs4)) ///
	ciopts(fcolor(gs13) lcolor(gs13) ) ///
	graphregion(color(white)) name(fp_weighted) ///
	xlabel(1(1)6, labsize(medlarge) format(%3.0fc)) xtitle("Price ($/Mcf)", size(large) margin(0 0 -2 2)) ///
	ytitle("Fitted Emission Rate", size(large) margin(0 2 0 0)) ///
	ylabel( `adjust_ax' "0" `=`adjust_ax'+.005' ".005"  `=`adjust_ax'+.01' ".01"  `=`adjust_ax'+.015' ".015" `=`adjust_ax'+.02' ".02", angle(vert) labsize(medlarge) glcolor(gs14) glwidth(vvthin) ) ///
	yscale(range(`=`adjust_ax'-0.003' `=`adjust_ax'+0.0201')) ///
	title("Weighting by Average Gas Production", size(large)) 

qui sum p_spot
graph addplot function y = `adjust_orig' + `original', lcolor(ebblue) lpattern("_##") lwidth(thick) range(`r(min)' `r(max)') ///
	ylabel( `adjust_ax' "0" `=`adjust_ax'+.005' ".005"  `=`adjust_ax'+.01' ".01"  `=`adjust_ax'+.015' ".015" `=`adjust_ax'+.02' ".02" `=`adjust_ax'+.025' ".025", angle(vert) labsize(medlarge) glcolor(gs14) glwidth(vvthin) ) ///
	yscale(range(`=`adjust_ax'-0.003' `=`adjust_ax'+0.0251')) ///
	saving(Graphs/fp_weighted, replace)

// Scatter showing facility size
qui sum p_spot
scalar min_p = `=r(min)'
scalar max_p = `=r(max)'
collapse (mean) rate gas p_spot, by(facility_id)
replace rate = .025 if rate > .025

twoway scatter rate p_spot [w=gas], msymbol(circle_hollow) mlwidth(vthin) mlcolor(gs8)

graph addplot function y = - `adjust' + `original', lcolor(ebblue) lpattern("_##") lwidth(thick) range(`=min_p' `=max_p') ///
	graphregion(color(white)) ///
	xlabel(1(1)6, labsize(medlarge) format(%3.0fc)) xtitle("Price ($/Mcf)", size(large) margin(0 0 -2 2)) ///
	ytitle("Emission Rate", size(large) margin(-2 3 0 0)) ///
	ylabel(0 "0" 0.005 ".005" 0.01 ".01" 0.015 ".015" 0.02 ".02" 0.025 ".025", angle(vert) labsize(medlarge) glcolor(gs14) glwidth(vvthin) ) /// yscale(range(-.001 .021))
	yscale(range(-.003 .0251)) ///
	legend(off) ///
	title("Average Facility Emission Rates by Size", size(large)) ///
	saving(Graphs/fp_size_scatter, replace)
	
graph combine Graphs/fp_weighted.gph Graphs/fp_size_scatter.gph, imargin(2 2 2 2) iscale(1) ///
	graphregion(color(white)) ///
	xsize(10) ysize(4)
	
graph export "$DATA/Results/p8_weights_size.pdf", replace



// 8.2 // Estimate multivariate FP model with interaction for price and gas

use "$DATA/ghgrp_di_intermediate", clear
drop if trim_5 == 1

gen p_spot_gas = p_spot*mean_gas

// First estimate linear model, then start with FP for coefficient with lower p value, then use estimated FP transformations for that one as controls in second one, then iterate until FP powers don't change
eststo: xtreg `dependent' p_spot p_spot_gas `controls' i.region_year `conditions' `weights', fe vce(cluster facility) // p_spot has lower p-value

fp <p_spot>, replace dimension(2) powers(-2 -2): xtreg `dependent' <p_spot> p_spot_gas `controls' `FE' `conditions', fe vce(cluster facility) // FP coefficients are -2 -2
gen p_fp1 = p_spot^(-2)
gen p_fp2 = p_spot^(-2)*log(p_spot)

// Now estimate FP for p_spot_gas controling for p_spot FP
fp <p_spot_gas>, replace dimension(2) powers(-.5 -.5): xtreg `dependent' <p_spot_gas> p_fp1 p_fp2 `controls' `FE'  `conditions', fe vce(cluster facility) // FP coefficients are -.5 -.5
gen p_gas_fp1 = p_spot_gas^(-.5)
gen p_gas_fp2 = p_spot_gas^(-.5)*log(p_spot_gas)

// Now return to p_spot
fp <p_spot>, replace dimension(2) powers(-2 -1): xtreg `dependent' <p_spot> p_gas_fp1 p_gas_fp2 `controls' `FE'  `conditions', fe vce(cluster facility) // FP coefficients are -2 -1
replace p_fp2 = p_spot^(-1)

// Now back to p_spot_gas
fp <p_spot_gas>, replace dimension(2) powers(-.5 -.5): xtreg `dependent' <p_spot_gas> p_fp1 p_fp2 `controls' `FE'  `conditions', fe vce(cluster facility) // FP coefficients are unchanged at -.5 -.5, so process is complete
mat B = e(b)
qui sum p_spot_gas if trim_5 != 1
local avg_p_gas `=r(mean)'
local rate_midpoint_R = B[1,`=colsof(B)'] + B[1,1]*(`avg_p_gas'^(-.5)) + B[1,2]*(log(`avg_p_gas')*`avg_p_gas'^(-.5))
local adjust_ax `=`rate_midpoint_R' - `avg_rate''
di `rate_midpoint_R'
di `adjust_ax'

fp plot if p_spot_gas > 20000000, residuals(none) plotopts(msize(small) mcolor(green)) lineopts(lwidth(medthick) lcolor(gs4)) ///
	ciopts(fcolor(gs13) lcolor(gs13) ) ///
	graphregion(color(white)) name(fp_interact_2) ///
	xlabel(0 "0" 500000000 "500" 1000000000 "1000" 1500000000 "1500" 2000000000 "2000" 2500000000 "2500" 3000000000 "3000", labsize(medlarge) format(%3.0fc)) xtitle("Price*Gas (millions)", size(large) margin(0 0 -2 2)) ///
	ytitle("Fitted Emission Rate", size(large) margin(0 2 0 0)) ///
	ylabel(`=`adjust_ax'-.02' "-.02"  `=`adjust_ax'-.01' "-.01" `adjust_ax' "0" `=`adjust_ax'+.01' ".01" `=`adjust_ax'+.02' ".02" `=`adjust_ax'+.03' ".03"  , angle(vert) labsize(medlarge) nogrid glcolor(white) glwidth(vvthin) ) ///
	ymlabel(`=`adjust_ax'-.02'(.005)`=`adjust_ax'+.03', nolabel tlcolor(white) grid glcolor(gs14) glwidth(vvthin)) ///
	yscale(range(`=`adjust_ax'-.0201' `=`adjust_ax'+0.032')) ///
	saving(Graphs/fp_interact_2, replace)

// One more return to p_spot for the graph and table output
eststo: fp <p_spot>, replace dimension(2) powers(-2 -1): xtreg `dependent' <p_spot> p_gas_fp1 p_gas_fp2 `controls' `FE'  `conditions', fe vce(cluster facility) // FP coefficients are -2 -1
mat B = e(b)
mat FP = e(fp_fp)
local rate_midpoint_R = B[1,`=colsof(B)'] + B[1,1]*(`avg_price'^FP[1,1]) + B[1,2]*(`avg_price'^FP[1,2])
local adjust_ax `=`rate_midpoint_R' - `avg_rate''
local adjust_orig `=`rate_midpoint_R' - `rate_midpoint''

fp plot, residuals(none) plotopts(msize(small) mcolor(green)) lineopts(lwidth(medthick) lcolor(gs4)) ///
	ciopts(fcolor(gs13) lcolor(gs13) ) ///
	graphregion(color(white)) name(fp_interact_1) ///
	xlabel(1(1)6, labsize(medlarge) format(%3.0fc)) xtitle("Price ($/Mcf)", size(large) margin(0 0 -2 2)) ///
	ytitle("Fitted Emission Rate", size(large) margin(0 2 0 0)) ///
	ylabel(, angle(vert) labsize(medlarge) glcolor(gs14) glwidth(vvthin) )

qui sum p_spot
graph addplot function y = `adjust_orig' + `original', lcolor(ebblue) lpattern("_##") lwidth(thick) range(`r(min)' `r(max)') ///
	ylabel(`=`adjust_ax'-.02' "-.02"  `=`adjust_ax'-.01' "-.01" `adjust_ax' "0" `=`adjust_ax'+.01' ".01" `=`adjust_ax'+.02' ".02" `=`adjust_ax'+.03' ".03"  , angle(vert) labsize(medlarge) nogrid glcolor(white) glwidth(vvthin) ) ///
	ymlabel(`=`adjust_ax'-.02'(.005)`=`adjust_ax'+.03', nolabel tlcolor(white) grid glcolor(gs14) glwidth(vvthin)) ///
	yscale(range(`=`adjust_ax'-0.0201' `=`adjust_ax'+0.032')) ///
	saving(Graphs/fp_interact_1, replace)

graph combine Graphs/fp_interact_1.gph Graphs/fp_interact_2.gph, imargin(2 4 2 2) iscale(1) ///
	graphregion(color(white)) ///
	xsize(10) ysize(4)
graph export "$DATA/Results/p8_size.pdf", replace

esttab using "$DATA/Results/p8_size.tex", b(4) se(4) replace star(* 0.10 ** 0.05 *** 0.01) keep(p_*) booktabs scalars(r2_w)
a
*/


